%% Load the data - Select the File: 'ViscosityProfiles'
close all; clear all; clc;
Dir = {};
Dir{1} = '.\ViscosityProfiles\0pt5PEO.mat'
Dir{2} = '.\ViscosityProfiles\0pt82PEO.mat'

first = 4; last = 8;    % Middle third of channel
t = [0 10 20 30 40 50 60];
gradCalc_Start = 2;     % 10 min
gradCalc_fin = 7;       % 60 min
t_10x_5 = [5 15 25];    % time points to be filled in
x = []; Viscosity_profile = {};Grad_Eta_label = [];
for i = 1:length(Dir)
    load(Dir{i}); 
    xlim_max = 1000; x(i,:) = xPos + (500 - median(xPos));
    Viscosity_profile{i} = cat(1,mu{1,:});
    FitViscSlope = [];FitVisc = [];
    for j = 1:length(mu)
        Mu_Exp = []; p = [];
        Mu_Exp = flip(mu{1,j}(first:last))
        p = polyfit(xPos(first:last),Mu_Exp,1);
        FitViscSlope(j,1) = p(1); FitViscSlope(j,2) = p(2);
        FitVisc = polyval(FitViscSlope(j,:),xPos(first:last));
    end
    Grad_Eta_label(i,:) = mean(FitViscSlope(4:10,1)); % How gradients are defined in Fig 1
    Grad_Eta(i,:) = FitViscSlope(gradCalc_Start:gradCalc_fin,1);
    eta_0(i,:) = Viscosity_profile{i}(gradCalc_Start:gradCalc_fin,6); % Find center spatial viscosity (\eta_0)

    grad_eta_eta0(i,:) = Grad_Eta(i,:)./eta_0(i,:);
    p = []; p = polyfit(t(gradCalc_Start:gradCalc_fin),grad_eta_eta0(i,:),2);
    p_fit(i,:) = p(3)+p(2).*t(gradCalc_Start:gradCalc_fin)+p(1).*t(gradCalc_Start:gradCalc_fin).^2;

    p_10x_5(i,:) = p(3)+p(2).*t_10x_5+p(1).*t_10x_5.^2; % Calculate grad eta/eta_0 for x = 5, 15, 25
    grad_eta_eta0_5minInterval(i,:) = [p_10x_5(i,1) grad_eta_eta0(i,1) p_10x_5(i,2) grad_eta_eta0(i,2) p_10x_5(i,3)];
end